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Abstract. - The conventional equations of Brownian motion can be derived from the first 
principles to order A 2 = m/M, where m and M are the masses of a bath molecule and a 
Brownian particle respectively. We discuss the extension to order A 4 using a perturbation 
analysis of the Kramers-Moyal expansion. For the momentum distribution such method yields 
an equation whose stationary solution is inconsistent with Boltzmann-Gibbs statistics. This 
property originates entirely from non-Markovian corrections which are negligible in lowest order 
but contribute to order A 4 . 



Dynamical justification of equilibrium statistical mechanics is a long-standing problem 
which can be traced to Einstein's criticism of the statistical definition of probability of 
macrostates [1] . The renewed interest stems primarily from recent development in the theory 
of nonextensive systems for which the validity of classical Boltzmann-Gibbs (BG) statistics 
is not obvious. A number of dynamical models generating non-canonical distributions have 
been suggested lately [2] , but for a truly conservative Hamiltonian system deviations from the 
BG statistics have never been dynamically justified. On the other hand, the BG statistics 
itself has been recovered directly from dynamics only for a few simplified models. Even for a 
dilute gas it has been done only in the limit of pair collisions. The ^/-theorem can be readily 
proved when dynamics is described by the master equation with transition rates obeying the 
detailed balance condition [3] . The latter can be deduced from reversibility of the microscopic 
processes involved. However, the master equation is usually a Markovian approximation of 
an originally non-Markovian exact equation for a targeted variable. The detailed balance 
arguments do not necessarily apply in this will be discussed at the end of the Letter. 

A phenomenon where relaxation to equilibrium may be traced, under reasonable assump- 
tions, from the first principles is Brownian motion of a heavy particle (B particle) of mass 
M in an infinite bath of light molecules of mass m <C M . Using an appropriate projection 
operator technique one can get to order A 2 = m/M the Langevin equation (LE) and the 
corresponding Fokker-Planck equation which describe exponential relaxation of the particle's 
momentum to the Maxwellian distribution. However, a generalization to higher orders in A is 
not trivial and has been studied in detail only for the very special Rayleigh model where the 
B-particle interacts with ideal gas molecules via the hard-wall potential and the molecules do 
not interact with each other at all [3,4]. 
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In this Letter we discuss Brownian motion to order A in the general case of continuous 
interaction. The essential difference with the Rayleigh model is that the Langevin equation 
in general is non-Markovian. However, one may perform systematic expansion around the 
Markovian limit converting the LE into a local equation with appropriate corrections. To 
order A 2 these non-Markovian corrections are negligible, but must be retained to order A 4 . 
The main result of the Letter is that non-Markovian corrections generate additional terms in 
the A 4 -order Fokker-Planck equation which make it inconsistent with the BG statistics. 

Let us first outline the elements of the dynamical theory of Brownian motion [5,6]. The 
basic assumption is that not too far from equilibrium the relation of the particle momentum P 
and typical momentum of a bath molecule p m is given by the equipartition theorem, p m /P ~ 
A = -y/m/M, and therefore the scaled momentum of the B-particle p = XP is of order A . It 
is assumed that writing the equation of motion for the scaled momentum p rather than for 
the true momentum P enables one to extract the overall dependence of the problem on the 
small parameter A. Then starting with the exact equation of motion for the B-particle and 
using an appropriate projection operator and perturbation techniques one can obtain to order 
A 2 the non-Markovian LE for the particle's scaled momentum 

^ = -A 2 £ drM (r)p(t - r) + \F{t). (1) 

Here the "random" force F(t) depends on the coordinates and momenta of bath molecules 
and fluctuates rapidly with mean zero. It is assumed that F may be expanded in powers of A 

F{t) = F (t) + XFi (t) + X 2 F 2 (t) + --- (2) 

where the first term Fo(t) has a meaning of the force exerted by the bath molecules on the 
B-particle fixed in space. In other words, F (t) is a fluctuating pressure on the infinitely 
heavy particle. Only this term should be taken into account in the A 2 -order equation, and the 
fluctuation-dissipation theorem has the form 

M (t) = (/3/m)(F (0)F (t)), (3) 

where f3 is the inverse temperature of the bath, and the average is taken with respect to 
canonical distribution of the bath in the field of the B-particlc held fixed. Next, noting that 
p(t — t) = p(t) — J. dt'pit') = p(t) + 0(A) one finds that non-Markovian effects are not 
important at order A 2 , and therefore the LE can be written in the conventional local form 

^ = -A 2 7oP (t) + AFoW, (4) 

with 70 = J °° dtMo(t). Integrating this equation one finds that {p 2 (t)) relaxes exponentially 
to the equilibrium value m/ (3 in agreement with BG statistics. 

Let us also remind the reader how the Fokker-Planck equation (FPE) for the distribution 
function f(p,t), corresponding to the LE may be derived. Starting with a generic mas- 
ter equation, which is valid for an arbitrary Markovian process, one can follow the familiar 
procedure [3] to transform it into the Kramers-Moyal expansion 



df( P ,t) 



dt n\ V dp 

n— 1 



h(-lY{ a ^ f ^} (5) 



where the coefficients are 



a n (p) = ]im- (\p(t + r)-p(t)} n }. (6) 

T—>0 T 
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These coefficients may be evaluated integrating the LE (0J for small r, 

ft + T 

p(t + T)-p(t) = -\ 2 loP (t)T + \ dt'F (t'). (7) 



This formula corresponds to a coarse-grained description with a time resolution much shorter 
than the relaxation time of the particle's momentum r p and much longer than the correlation 
time for the random force r c . One can not require a finer time resolution because the LE 
itself in the form |0J is valid only on a time scale t » t c . Two remarks are in order here. 
First, the route we follow is not a good one when r c does not exist, for instance when C(t) = 
(Fo(O)Fo(i)} ~ t~ a ' . We ignore the slow hydrodynamic modes of the bath responsible for 
long-time tails of correlation functions. Second, the coarse-grained description implies that 
the moments ([p(t + r) — p(t)] n ) in Eq.@ must be calculated in the limit r 3> t c . Only after 
that must the operation lim T ^o h(- • ■) be taken. 

In general the "random" force in the LE is neither delta-correlated nor Gaussian, and 
therefore the Kramers-Moyal expansion is not truncated. However, only the two first terms 
involve contributions of order A 2 . Then to order A 2 the expansion turns into the familiar 
second-order Fokker-Planck equation 

dJ ^ 1 = X 2 D 2 f( P ,t), D 2 = ai ^p + a 2 ^. (8) 

Since the pressure F (t) is a stationary process, the coefficients have the form 

Oi=7b = -/ dt(F (0)F (t)), a 2 = i / dt x [ dt 2 (F (*i)F (<2)>. (9) 

Recall that the double integral in the last expression must be taken in the limit r 3> t c . The 
subsequent limit r — > is not necessary because the integral grows linearly with time for 
t ^> t c . Indeed, taking into account that (F (ti)Fo(t 2 )) = C{\t\ — t 2 \) one obtains 

pT pT pT pt\ pT 

dtJ dt 2 (F (t 1 )F {t 2 ))=2 dt x dt 2 C(h-t 2 ) = 2 dti(T-ti)C(ti). (10) 
Jo Jo Jo Jo Jo 

For r 3> t c this integral equals 2r §^ dt\C(t\). With this result, one finds the relation 

Oi = — a 2 . (11) 

m 

This means that the stationary solution of the equation (JSJ, f(p) — cexp^— ^-p 2 ^, is 
Maxwellian. 

Our aim is to generalize the above approach to order A 4 . (The A 3 -order damping force 
vanishes due to symmetry for a homogeneous bath.) The A 4 -order LE has the form 

^^ = -A 2 / dr M x {T)p(t - t) - A 4 / d T M 2 (T)p 3 (t~T) + XF(t). (12) 
dt J J 

A microscopic derivation of this equation was recently discussed in detail in [6] . It differs from 
the A 2 -order LE Q not only by the presence of the nonlinear damping term involving p 3 , but 
also by additional corrections of order A 2 to the memory kernel for linear damping 

Mi(t) = M (t) + A 2 SM(t). (13) 
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Another difference is that the "random" force in Eq. i|12|) can not be approximated by the 
pressure term Fo(t), as in the linear LE QJ, but must involve the corrections of higher orders 
F\ and F 2 , see Eq. (J2J) - The correction to the linear damping SM and the nonlinear damping 
kernel M 2 (t) can be expressed in terms of correlation functions involving F$, F\ and F 2 [6], 
but their explicit forms are not important for our purpose here. 

Our first step is to find a Markovian approximation for the non-Markovian equation l)12|l 
using the same trick as for the linear LE, i.e. writing (j>{t — r) = cf>(t) — J? dt'<j>(t') for 
4>(t) = p n {t). Since -^.p n (t) ~ pit) = 0(A) according to 1|12[1 . the non-Markovian correction 
for the nonlinear damping is of order A 5 and should be neglected, 



—A / dr M 2 (r)p 6 (t - t) -» -AVW / dr M 2 (r) (14) 
Jo Jo 

However, for the linear damping term the non-Markovian corrections are of order A 3 and A 4 
and must be retained. Indeed, let us write the linear damping term as a local expression plus 
a correction A(t), 



-X I drAfi (r)p(i - r) = -X 2 p(t) drMi(r) + A(t). (15) 
Jo Jo 

The correction 

A(£) = A 2 / drMi(r) / fft'p(t') (16) 

JO Jt-T 

can be evaluated estimating p(£) from the LE (|12|l . For i 3> r c one obtains to order A 4 

A(i) = A 3 / dr M (r) / dt'F {t')-X A p{t) I dtM (t) / dtM (t)t. (17) 
Jo Jt-r Jo Jo 

We shall write this expression in the form 

A(*) = A 3 ^(i)-AVp(*)- (18) 
Here the first term involves a rapidly changing function 

F*(t) = 2- [ dr C(t) f dt'F (t') (19) 

m JO Jt-T 

and may be considered as an additional contribution to the random force. Note that F*(t) is 
approximately stationary on the time scale t ^> t c . The second term in Ea. H18(l with 

7 * = (m) j( dt °^ Jo dtC ^ t ^ 

can be interpreted as an additional contribution to the linear damping. 
As a result, the LE (|12f> can be written in local form as follows 

^ = -X 2 llP (t) X^ 2P \t) + A£(t). (21) 

Here is the "random" force with a A 2 -order non-Markovian correction 

£(t) = F(t) + X 2 F*{t) = F (t) + \Fi(t) + X 2 F 2 (t) + X 2 F*{t). (22) 



A. V. Plyukhin: Does a Brownian particle equilibrate? 



5 



The linear damping coefficient in l|21|) takes the form 

7i = 7 o + A 2 <5 7 + AV, (23) 

with 5j = Jg 00 dt5M(t), and the nonlinear damping coefficient is 72 = / °° dtM 2 (t), 
Integrating the LE (|2*l"1) for r c <C r <C t p one gets 

p(t + r)~p(t) = -[X 2 ll p(t) + X 4 1 2P 3 (t)] r + \j 'dt'Zit'). (24) 

Then using Eq. (jfJJl one finds the coefficients a n in the Kramers-Moyal expansion J5J , 

Q'i = — A 2 7i p - A 4 72 p 3 — - A 2 7o p - A 4 7* p — A 4 (£7 p + 72 p 3 ) (25) 
and, recalling that £(£) is stationary for t 3> r c , 



« 1 



= X n lim- [ dh f dt 2 ... f dfcn<£(*i)£(*a)-£(*»)> ( 26 ) 
r Jo Jo Jo 



for n = 2, 3, 4. (The terms with n > 4 do not contribute at order A .) In the above expression 
the integrals must be taken in the limit t 3> r c . 

According to l|26|l . to get 02 to order A 4 one needs the correlation (£,(ti)!;(t 2 )) to order A 2 , 

<£(*i)£(*a)> = <Fo(ti)F (t 2 )> + A 2 <Fo(ti)F*(t 2 )) + A^Fq^F*^)) 
+A 2 (F 1 (t 1 )F 1 (i 2 )) + A 2 (F (ii)F 2 (t 2 )) + A 2 <F (t2)F 2 (ti)>. 

Similarly, a 3 and 04 require the correlations to order A 1 and A , respectively: 

(£(ti)£(t2)£(t 3 )) = A{F (*i)Fo(i2)Fi(t 3 )} + A(^o(ti)Fi(i 2 )Fo(i 3 )) + A{Fi(ii)Fo(t 2 )Fo(i 3 )), 
(^1)^2)^3)^4)) - (^o(ii)^ofe)^ofe)^o(i4)). 

Using these formulas and collecting the terms of the same order, one obtains the A 4 -order 
FPE in the form 



df(p,t) 
dt 



= [x 2 D 2 + X A D* 2 +X i D i )f{p,t). (27) 



Here the differential operator D 2 is the same as in the A 2 -order FPE 1(8} . The operator D 2 
originates from non-Markovian corrections and has the same structure as D 2 , 



D * =b 4 P P + b *$ < 28 > 



with b\ = 7* given by Ea. (j2()|) and 



l -f f dt 1 dt 2 {(F*(t 1 )F (t 2 )) + (F*(t 2 )F (h))} = ^ f [ dt 1 dt 2 {F*(t 2 )F (h)). (29) 



h , , , 

2t Jo Jo 

The operator D4 in Eq.J2Zjl is a differential operator of order four. It absorbs the terms 
of order A 4 which are expressed in terms of time correlations (FqF 2 ), (F\Fi), (FqFoFi), and 
(FqFqFqFo) . Since (Fi(t)) = 0, the first three of these correlations are actually cumulants, 
which we denote by ((• ■ ■)), 

(FiFj) ee (Ej,)(Fj) + ((FiFj)) = ((FiFj)), 
(FiFjFk) ee (F) (Fj)(F k ) + (Fi)({FjFk)} + <F,-)<(FF fc }) + (F fe >«i^» + ((FF^)) = ((FF^)). 



6 



EUROPHYSICS LETTERS 



The cumulant expansion of the correlation (FqFoFqFo}, which determines a^, involves the 
product of cumulants ((Fo(ti) Fofo))) ((-^0(^3) -^0(^4))} and two similar terms. However, since a 
product of cumulants depends only on two time differences, the corresponding contributions 
to the four-dimensional time integral in the expression for «4 vanish in the limit r — > 0. 
Therefore only the cumulant ((^0(^1)^0(^2)^0(^3)^0(^4))) contributes the expression for 04. 

The fact that D4 is linear in cumulants of Fi means that D4 is a linear functional of 
cumulants for the density of bath particles N(z,t) = 8(z — Zi(t)), where zi denotes the 
coordinate-momentum pair (xi,pi) of a bath particle. In turn, one can observe that cu- 
mulants ((N(zi,ti)N(z2,t2).--N(zk,tk))) of any order k depend linearly on the concentra- 
tion of bath molecules n. For instance, in the expression for the product N(zi, ti)N(z2, £2) 
one can write the double sum as J2i j = Si^ +Si=j which gives {N(zi,ti)N(z2,t 2 )) = 
(N(zx,ti))(N(z2,t2)) + J2i(d( z i ~ z i{ti))${ z 2 — Zi(t2))). Here the second term on the right 
side is by definition the cumulant {{N(zi,ti)N(z2,t2))) which is obviously linear in n. 

It follows from the above discussion that the operator depends linearly on the concen- 
tration of bath molecules n. In contrast, the first and second terms in D\ are proportional 
respectively to 7* and the correlation (Fo(t\)F*(t2)), which are both quadratic in n. The 
Maxwellian distribution /m(p) ~ exp (— /3p 2 /2m) depends neither on A nor on n. Therefore, 
if /m(p) is a stationary solution of the FPE (|27|) . it must satisfy each term separately: 

D 2 fu(p)=0, D* 2 f M (p)=0, D i f M (p) = 0. (30) 

The first relation, as was discussed at the beginning of the Letter, is satisfied. The validity of 
the relation D^/m — is not our concern here. One can show that -D4 has the same structure 
as for the Rayleigh model [3,4] and the validity of the relation D4/M = can be explicitly 
checked for the extended Rayleigh model [6]. The central result of this Letter is that the 
second equation D^Jm — is not satisfied. Indeed, this equation implies that b\ — —62, 
similar to the condition Ullfl for Z?2- Instead, the relation with the additional factor 1/2 
actually holds: 

6l = A 6 , (31) 

Explicitly this relation has the form 

f-OO pOC 1 pT pT 

/ dtC(t) dtC{t)t = — / dt 1 dt 2 C*(t 1 ,t 2 ), (32) 
Jo Jo 2r J Q J 

where C{t) = (F o {0)F Q {t)) and C*(h,t 2 ) = (fi/m)- 1 (F Q {t{)F* (t 2 )) . According to (JTHJ), th e 
latter correlation is given by 

C*(h,t 2 )= r dt'C(t') r dt"C{\t 1 -t"\). (33) 
Jo Jt 2 -f 

Note that Eq. (|32J) is an asymptotic relation: the double integral of C* must be taken in the 
limit t ^ t c . The validity of the relation (|32|l may be verified directly, for instance, for 
C(t) ~ e~*/ Tc or for a Gaussian correlation. A general proof of (|32|l is rather lengthy and will 
be presented elsewhere. 

Since D^fnip) 7^ 0, the stationary solution of the FPE (|27|l is not Maxwellian and in 
general can not be analyzed without identifying an explicit form of the operator D4 which is 
beyond our scope here. Suppose however that the concentration of bath molecules n is in a 
sense large and -D4 ~ n can be neglected compared to D\ ~ n 2 . The A 4 -order FPE with £> 4 
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dropped is an equation of second order which has the same structure as the A 2 -order equation 
iJHJ but with modified coefficients a* — cii + X 2 bi. Its stationary solution exp ^rp 2 ^j is 

the Maxwellian distribution exp (^~-§mP 2 ^j but with the inverse temperature (3* — ma\ja* 2 
slightly lower than that for the bath, f3. Taking into account relations and i)31|l . one 
obtains /3* = (3 (f — A 2 b\/a\), or 

0* = /?jl-A 2 £ J dtC(t)t\. (34) 

This result may be interpreted alternatively as a renormalization of the mass of the B particle 
M — > M* — {p/p*)M. 

In conclusion, in this Letter we discussed the derivation of the Fokker-Planck equation 
corresponding to the microscopic Langevin equation of order A 4 . The procedure includes the 
expansion of the originally non-Markovian Langevin equation around the Markovian limit and 
evaluation of the coefficients in the Kramers-Moyal expansion of the corresponding master 
equation. Let us emphasize that the method we used here is not just a truncation of the 
Kramers-Moyal expansion, which is known to yield an incorrect description. Rather, we 
transformed the Kramers-Moyal expansion into the expansion in powers of small parameter 
A and applied a perturbation technique to obtain an equation to given order in A. This 
may be seen as a generalization of van Kampen's 1/il-expansion [3,4]. It is found that to 
order A 4 the method leads to the equation which has a non-Maxwellian stationary solution. 
This feature originates entirely from non-Markovian corrections collected in the operator D 2 , 
Eq.J2Hl- Since the method yields correct results for the Markovian Rayleigh model [4] at 
least to order A 4 , one may conclude that either the deviation from Maxwellian distribution for 
conservative Brownian motion is a small but real effect, or the procedure adopted here does 
not treat non-Markovian effects properly. Note that we assumed but did not prove that each 
term in A-expansions is bounded for all t. If the assumption is correct, the results of the Letter 
would mean that the conventional time-reversal arguments [3] , leading to the detailed balance 
and, eventually, to the H-theorem, do not apply to Brownian motion beyond the lowest order 
in A. This is perhaps not such a surprise recalling that the right part of a Langevin equation 
is not a true force on the particle but rather an approximation to given order in A. Whether 
or not this truncated force corresponds to any effective Hamiltonian is not known a priori, 
and therefore the invariance for time reversal can not be invoked. On the other hand, an exact 
closed equation for a targeted variable is inevitably non-Markovian, and the standard detailed 
balance arguments can not be applied either. As a referee noted, there is a similarity with the 
Onsager relations which follow from microscopic reversibility but become only approximately 
valid when higher order corrections are included in a fast variables elimination scheme [7]. 
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